OPEN 3 ACCESS Freely available online 



•0-PLOS I o-^E 



Scaling Up Stomatal Conductance from Leaf to Canopy (8) 
Using a Dual-Leaf Model for Estimating Crop crose^rk 
Evapotranspiration 

Risheng Ding\ Shaozhong Kang^% Taisheng Du\ Xinmei Hao\ Yanqun Zhang^ 

1 Center for Agricultural Water Research in China, China Agricultural University, Beijing, China, 2 National Center of Efficient Irrigation Engineering and Technology 
Research-Beijing, China Institute of Water Resources and Hydropower Research, Beijing, China 

Abstract 

The dual-source Shuttleworth-Wallace model has been widely used to estimate and partition crop evapotranspiration {).ET). 
Canopy stomatal conductance (GjJ, an essential parameter of the model, is often calculated by scaling up leaf stomatal 
conductance, considering the canopy as one single leaf in a so-called "big-leaf" model. However, G^^ can be overestimated 
or underestimated depending on leaf area index level in the big-leaf model, due to a non-linear stomatal response to light. 
A dual-leaf model, scaling up Gj^ from leaf to canopy, was developed in this study. The non-linear stomata-light relationship 
was incorporated by dividing the canopy into sunlit and shaded fractions and calculating each fraction separately according 
to absorbed irradiances. The model includes: (1) the absorbed irradiance, determined by separately integrating the sunlit 
and shaded leaves with consideration of both beam and diffuse radiation; (2) leaf area for the sunlit and shaded fractions; 
and (3) a leaf conductance model that accounts for the response of stomata to PAR, vapor pressure deficit and available soil 
water. In contrast to the significant errors of C^c in the big-leaf model, the predicted G^c using the dual-leaf model had a high 
degree of data-model agreement; the slope of the linear regression between daytime predictions and measurements was 
1.01 (R^ = 0.98), with RMSEofO.6120 mm s for four clear-sky days in different growth stages. The estimates of half-hourly 
aET using the dual-source dual-leaf model (DSDL) agreed well with measurements and the error was within 5% during two 
growing seasons of maize with differing hydrometeorological and management strategies. Moreover, the estimates of soil 
evaporation using the DSDL model closely matched actual measurements. Our results indicate that the DSDL model can 
produce more accurate estimation of G^c and lET, compared to the big-leaf model, and thus is an effective alternative 
approach for estimating and partitioning lET. 
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Introduction 

Accurate estimation of evapotranspiration [aET) is important in 
understanding terrestrial hydrological cycles because XET is the 
largest component in the terrestrial water balance after precipi- 
tation [1]. In agricultural production, improved estimation of crop 
lET is also needed to develop precise irrigation scheduling and 
enhance water use efficiency, as soil water depletion is mostly 
determined by the rate of aET [2,3,4]. However, direct 
measurement of XET is often difficult, costly and not available in 
many regions [5,6]. Therefore, mathematical models are needed 
to estimate lET using readily measurable meteorological and 
environmental variables. 

Vegetation transpiration (T,) and soU evaporation (£j), which 
are controlled by different biotic and physical processes, are the 
two major components of XET. Transpiration is strongly linked to 
crop productivity since it occurs concurrendy with photosynthetic 
gas exchange [7]. Quantifying T,- is also critical to accurately 
predict the response of crop functioning and physiology to 



changing climate [8]. Because the two separate processes occur 
simultaneously, there is no simple way to distinguish between them 
[9,10]. 

Several models have been developed to calculate XET and 
separately estimate soil evaporation and transpiration [11,12,13]. 
Shutdeworth and Wallace [14] described a dual-source model 
with a resistance-energy combination, which could separately 
predict T,- and £j, and is also sufficientiy simple [8]. This model 
has been widely used and can also be used to gain an 
understanding of the interaction of biophysical and hydrological 
processes in the crop canopy [14,15]. Determination of different 
resistances or conductances (the reciprocal of resistance) is 
necessary for its practical application. Specifically, canopy 
stomatal conductance (G„) is often calculated by scaling up leaf 
stomatal conductance of the leaves acting in parallel while treating 
the canopy as one big-leaf, hereafter the big-leaf approach [8, 1 6] . 

The weakness of using the big-leaf approach is that the use of 
mean absorbed radiation can significantly overestimate G„, 
especially in dense canopies, because the light response of stomata 
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measurements taken in an irrigated maize field mulched with 
plastic film. 

Model Development for Evapotranspiration and 
Canopy Stomatal Conductance 

1 Dual-source Evapotranspiration Model 

The XET in the dual-source model was partitioned into two 
components, canopy transpiration (kT^ and soil evaporation [XE^ 
with a resistance network [14]. 



kET = XTc + kE, = cOcPMc + m,PM, 



(1) 



^ AA + {pCpVPD-Ara,A,)/{raa + ra,) 

A + y[l + r,c/{raa + rac)] ^ ' 



Figure 1. Response of canopy stomatal conductance (C^^) to 
incident PAR above the canopy (C?J at different leaf area 
indices [LA/). Gj^, and G;^2 are calculated by the big-leaf and dual-leaf 
canopy stomatal resistance models, respectively. 
doi:1 0.1 371/journal.pone.0095584.g001 

is non-linear [17,18]. Moreover, the overestimation of G„, by the 
big-leaf approach often occurred when using total leaf area index 
(LAT) to scale up stomatal conductance. To mitigate the 
overestimation of G,„ researchers have used the effective or sunUt 
LAI (LAIg) instead of total LAI [16,18]. But, the relationship 
between LAI^ and LAI is empirical and varies with the vegetation 
type and solar radiation angle. Therefore, models of Gsc should 
consider the non-linear response of stomata to irradiance as well as 
heterogeneous radiation profiles in the canopy. Beam and diffuse 
radiation penetrating the canopy must be considered separately 
due to differential attenuation in canopies, as should visible and 
near-infrared wavebands due to differential absorption by leaves 
[12,19]. A multilayers model could account for the mechanism of 
radiation penetration [20,21], but it is useless in practical 
applications [22]. Some studies indicate that radiation penetration 
within the canopy could be simplified by splitting the canopy into 
two fractions of leaves: sunlit and shaded [18,23,24]. Leuning et al. 
[25] used a multilayer model to show that photosynthesis from a 
canopy is closely approximated when calculated as the separate 
sums of sunlit and shaded fractions, weighted by their respective 
leaf area within the canopy. However, few studies have calculated 
using the approach of sunlit and shaded leaves (hereafter, the 
dual-leaf approach) in the dual-source S-W model for estimating 
XET. 

Soil surface resistance is another key parameter for 
partitioning ?.ET in the dual-source model because above-canopy 
XET\% contaminated by flux from the soil substrate as variations in 
leaf area affect soil exposure, soil evaporation and absorbed 
radiation [14,26]. To reduce soil evaporation in the field, the 
ground is often mulched with plastic film, and this technique is 
widely used in northwest China [10,27]. But the effect of ground 
mulching on soil evaporation has not been taken into account 
when param("t("rizing r,, in th(" dual-source model. In this study, a 
dual-leaf model of G„. was developed by scaling up leaf stomatal 
conductance using the corresponding absorbed irradiance for the 
sunlit and shaded leaves separately, which overcomes the 
Umitations of the big-leaf model. The dual-leaf model developed 
here was then incorporated into the dual-source model to estimate 
and partition lET. We evaluated the dual-source dual-leaf model 
(DSDL) by comparing the model estimations of kET with 



^ AA + [pCp VPD-Ar^jA-A,)] /(r,, + r,,) 



A + y[\-lrr,,/{raa + ra. 



1 

\+R,Ral[R,{Rc+Ra)\ 



\+R,Ra/[Rc{Rs + Ra) 



Rs={A+y)ras+yrs, 



Rc = {A + y)rac + yrsi 



Ra = {A + y)ra, 



(4) 



(5) 



(6) 



(7) 



(8) 



where PAf, and PM^ are the terms similar to those in Penman— 
Monteith model for canopy transpiration and soil evaporation, 
respectively, and Wc and cOj are the weighting factors for the crop 
canopy and soil components, respectively. X is the heat of water 
vaporization, p is air density, Cp is the specific heat of dry air at 
constant pressure, A is the slope of the saturation vapor pressure 
curve, y is the psychrometric constant, VPD is vapor pressure 
deficit, and A and A, are the total available energy and available 
energy for soil, respectively, r^c is canopy stomatal resistance, r,, is 
soil surface resistance, Tac is canopy boundary layer resistance, r„, is 
soil boundary layer resistance between soil and vegetative canopy, 
and Taa is aerodynamic resistance between canopy source and 
reference height, respectively. The calculation procedures of the 
other resistances except r^c and tss are given in Appendix SI. 



A = R„-G 



A, = R„ 



(9) 



(10) 
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7?^ = R„ exp ( - krLAI) 



(11) 



where R„ and are net radiation al)o\c tiic canopy and at the 
soil surface, respectively, and G is the soil heat flux. The canopy 
extinction coefficient of net radiation, Kji, is dependent on leaf 
orientation and solar zenith angle {Q [17]. 



2 Irradiance within Crop Canopy 

Incident PAR light above the canopy {Q^ was divided into 
diffuse {Q^j) and beam irradiance ((^j) through the fraction of 
diffuse radiation {fj). 



Qod=fdQo 



(18) 



Kr = 



cos(0 



(12) 



where Gj^ is 0.5 for a spherical leaf angle distribution. ^, the angle 
subtended by the sun at the center of the earth, is perpendicular to 
the surface of the earth and calculated as in Appendix S 1 . 

The two components, XT^ and XE^ were now calculated along 
with the VPD at the canopy source height (D„)- 



D„=VPD+[AA-{K+y)XET\raa/ pCp 



XT, 



A{A-A,)+pCpDo/ra. 
A + y(l+r,,/ra,) 



XEs- 



(13) 



(14) 



(15) 



The measured r^^ was obtained by inverting Eq. (14), with XT^ 
calculated by the known or measured XET and XE^. 



X Tr = XET — XE^ 



(16) 



Qob = {i-fd)Qo 



(19) 



The f,/ was cak:ulated from a simple model of atmospheric 
attenuation of radiation [17,28]. 



1-T,,"'" 



i+T„™«(i/./;,-i) 



(20) 



where t„ is the atmospheric transmittance, is the forward 
scattering coefficient of PAR in atmosphere, and is the optical 
air mass, which can be calculated as follows. 



P/Pq 
cosC 



(21) 



where P is loc al atmospheric pressure and Pq is atmospheric 
pressure at sea level. 

At a depth c in the canopy, three types of irradiance can be 
calculated: the total beam, Qf^^, (unintercepted beam plus down 
scattered beam), direct beam, j (unintercepted beam) and the 
diffuse flux, Qf^^ [17,29]. 

Qe,biiO = Qob{l-Pcb)V^Ki, exp(- v^KftiJ) (22) 



A{A-As)+ PgCpDolrg, - XT,{A + y) 
yXTc/rac 



(17) 



Qe,b(.0 = Qob{i-Pcbhb exp( - Kbi) (23) 
QeAi) = Qod{^- Pcd)V«Kd exp( - VaKdi) (24) 



PcA = 1 - exp(2p^(Ci/(l + Kb)) 



(25) 




Figure 2. Variation in the sunlit and shaded leaf area indices 

{LA/s/ and LAIs/,) versus different LAI. The ratios of LAI^i and LAI^t, to 
LAI (Asi and A^^i are also shown. 
doi:1 0.1 371/journal.pone.0095584.g002 



Ph = 



1 -i/o( 
1 + Va 



Pcd = 



2KdP^ 
Kj+l 



(26) 



(27) 



where a is absorptivity of leaves for irradiation, p^/, and p^j are 
canopy reflectance for beam and diffuse irradiance respectively 
with a randomly spherical leaf-angle distribution, p/, is canopy 
reflectance for beam irradiance with a horizontal leaf-angle 
distribution, and Kj is an extinction coefficient for difiuse radiation. 

The absorbed irradiance in a canopy height {Q/j consists of the 
total beam radiation {Q^^i,^ and the diffuse radiation (Q^^J). 



m)=Qe,bM)+Q(AO 



(28) 
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The total irradiance absorbed by the entire canopy (Q^) per unit 
ground area was determined by integrating over the total LAI. 



Qc= I Qmdi=Qob{\-p,b)[^-exv{-^/^KbLAI)] 

0 



(29) 



The irradiance absorbed by the sunlit fraction in a specific: 
canopy height {Q^^ij can be given as the sum of direct-beam (Qf^j), 
diffuse (Q^^rf) and scattered-beam components [QiJ. 



QeAi) = QeAO + QeAO+Qe.M) 



(30) 



The total irradiance absorbed (Q^) is the sum of the two parts, 
irradiance absorbed by the separate sunlit (Q^l) and shade fractions 
iOsh) of the canopy. Thus, Q^/, was calculated as the difference 
between and Q,;. 



Qsh = Qc — Qsl 



(33) 



3 Leaf and Canopy Stomatal Conductance 

The stomatal conductance is represented by g, for a single leaf 
and G„. for the entire canopy. 

3.1. Leaf stomatal conductance. The leaf can be 
calculated using the Jarvis-Stewart type multiple formulae 
[30,31,32]. 



QUO 

= Qob [(1 - Pcb)V«'Kb exp( - VaKb^) - aKb exp( - Kbi)] 



(31) 



The irradiance absorbed by the sunlit fraction in the entire 
canopy {Qjij was obtained by integrating j; over the total LAI. 



0 

LAI LAI 

= I Qej,{OkM)d^+ I QiAOkM)d^ (32) 



r Ssmax ^ i ^X: 



(34) 



where ^^m, is the maximum value of the leaf stomatal conductance 
and Fxi is the stress function of the specific environmental variables 
(.r;), O^Fx,^l. The original model used short-wa\'e radiation as the 
light variable. Here we have used the photosynthetically active 
radiation absorbed by canopy leaves {Q^ because stomatal 
aperture is determined by the received visible wavelength 
radiation, rather than short-wave radiation [17,18,33]. In addi- 
tion, we incorporated the environmental stress impact on g^ by 
VPD and available soil water as follows. 



500 

500 Qa+kQ 



(35) 



(36) 



I Qi,b(S)fiMiOd^ = QobOi[l - exp( - KbLAI)] (32a) 



1 - exp{-k,r6E) 
l-exp{-k^) 



(37) 



LAI 

I QeAiVeAOd^ 



(32b) 



y^Kd 



QodO--Pcd)—r \ [1- exp(-(VaKrf-|-KA)L^/)] 

\/a.Kd + Kb 



LAI 

I QtAOfUOdi 



0 

= Qob 



(32c) 



'•^^'Pcb) (1 - exp ( - yiKb -I- Kb) LAI) ) 



-a{\-s%.p{-2KbLAI))ll 



Op — On 



(38) 



where the Rq, ko and k„ are the stress coeflicients of Q^, VPD and 
extractable soil water in the root zone (0^), and 6, 9p and Ojy are 
the measured soil moisture, field capacity and wilting point in the 
root zone, respectively. 

3.2. Big-leaf model of canopy stomatal 
conductance. The canopy stomatal conductance in the big-leaf 
model (Gsci) is estimated by scaling up gs weighing by the effective 
LAI [LAI^ as if the canopy is a single big-leaf [16,18,32]. 



Gsc\ = =gs\LAIe 



(39) 



where g,i is the mean leaf stomatal conductance for the entire big- 
leaf and can be calculated by Eq. (34) based on the mean absorbed 
irradiance of the entire canopy using Eq. (29); LAI, is empirically 
equal to the actual LAI tor LAI<2, LAI/2 for Z4/>4, and 2.0 for 
others [16]. 
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Table 1. The paranneters used in the dual-source dual-leaf model. 



Symbol 


Description 


Value 


Units 


Sources 


b, 


Parameter in soil resistance model 


15.2 


sm-^ 


Fitted in this study 


b2 


Parameter in soil resistance model 


-5.8 




Fitted in tiiis study 


bs 


Parameter in soil resistance model 


88.7 




Fitted in this study 




Mean drag coefficient for the individual vegetative elements 


0.1 




Meyers and Paw [42] 


d, 


Characteristic leaf dimension 


0.068 


m 


Measured in this study 




Spherical leaf angle distribution 


0.5 




Campbell and Norman [17] 


Qsmax 


Maximum value of the leaf stomatal conductance 


7.5 


mm s"^ 


Fitted in this study 


k 


von Karman's constant 


0.41 


- 


Brutsaert [41] 


kg 


Stress coefficients of the 
photosynthetically active radiation 
in the stomatal 
conductance model 


150 


W m"^ 


Fitted in this study 




Stress coefficients of the vapor pressure deficit in the stomatal conductance 
model 


0.2 


kPa"' 


Fitted in this study 


k„ 


Stress coefficients of the extractable soil water in the root zone in the stomatal 
conductance model 


7 5 




Fitted in this study 


Zos 


Effective roughness length of the soil substrate 


0.01 


m 


Shuttleworth and Wallace [14] 


oc 


Absorptivity of leaves of irradiation 


0.8 




Monteith and Unsworth [12] 




Extinction coefficient for diffuse radiation 


0.7 




Campbell and Norman [17] 




Atmospheric transmlttance 


0.72 




Brutsaert [41] 


fa 


Forward scattering coefficient of PAR In the atmosphere 


0.43 




Brutsaert [41] 


frr, 


Fraction of ground mulched by plastic film 


0.5/0.6= 




Measured in this study 



"0.5/0.6 is the fraction of ground mulched by plastic film in 2009 and 2010, respectively. 
doi:l 0.1 371 /journal.pone.0095584.t001 



3.3. Dual-leaf model of canopy stomatal conductance. In 

the dual-leaf model, G„ [G„2) is calculated by summing the 
contributions of sunlit and shaded fractions, G^; and G^k, 
respectively, which are scaled up using the associated g, weighted 
by their respective fractions oi LAI [18,24]. 



LAL, 



f,M)d^='^^'^^^'^ (43) 



1 

rsc2 



-- Gsi + Gsh =g,iLAI,i+g,i,LAI,i, 



(40) 



where g^i and g,/, are tiie mean leaf stomatal conductance for sunlit 
and shaded leaves, respectively, and can be calculated by Eq. (34) 
based on the separate absorbed irradiance using Eqs. (32) and (33). 
LAI,i and LAI,/, are LAI for sunlit and shaded leaves in the entire 
canopy, respectively. 

Assuming that all leaves in a canopy are randomly distributed, 
the fraction of sunlit leaves {f(^,i) in a specific canopy depth declines 
exponentially with cumulative leaf area {£) [29,34]. 



fiAO= exp(-KhO 



Ld(z)dz 



(41) 



(42) 



where is an extinction coefficient for beam radiation, Lj is leaf 
area density, z is height above ground, and he is canopy height. 
LAIsi is therefore calculated by integrating ,/ for the entire 
canopy. 



LAL,, = LAI -LAL 



si 



(44) 



4 Soil Surface Resistance 

In this study, r„ was directly calculated with a function 
dependent on surface soil water content [35], accounting for the 
effect of plastic mulching on reduction of soil evaporation by 
introducing a term for fraction of plastic mulch, f„, [i.e. r„ is 
divided by the area of exposed substrate per unit ground area ( 1 — 



1 



(45) 



where Og is the average soil water content between 0-0. 1 m, 6, is 
the saturated water content of surface soil, and hj, A^, and 65 are 
the empirical coefficients. 
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Figure 3. Diurnal variations of (a) modeled irradiance, (b) leaf 
area index and (c) canopy stomatal conductance at leaf area 
index [LA/I = 5.0. Qo is the total irradiance above the canopy, is the 
irradiance absorbed by the entire canopy, and is separated into the 
irradiance absorbed by the sunlit leaves of the canopy (Qs/) and 
irradiance absorbed by the shaded leaves of the canopy (Qj/,). LAI;/ and 
M/j/, are the sunlit and shaded fractions of LAI, respectively. G^d and 
Gsc2 are the canopy stomatal conductance calculated by the big-leaf 
and dual-leaf models, respectively; 65^2 is separated into the sunlit and 
shaded canopy stomatal resistance, Gsi and Gsh, respectively. 
doi:1 0.1 371/journal.pone.0095584.g003 

Materials and Methods 

1 Experimental Arrangement 

The experiments were conducted in an irrigated maize field in 
2009 and 2010 at Shiyanghe Experimental Station for Water- 
saving in Agriculture and Ecology of the China Agricultural 
University in Gansu Province in northwest China (N 37°52', E 
102°50', Altitude 1581 m). Grain maize was sown on April 21 and 
May 2, and harvested on September 28 and September 26 in 2009 
and 2010, respectively. The ground was partly mulched with 
plastic films with a width of 1 00 cm, and there was bare soil of 65 
and 45 cm in width between two film rows in 2009 and 2010, 
respectively. Maize seed was sown in holes of 5.0 cm diameter 
under the plastic films, with a row spacing of 50 and 45 cm and a 
plant spacing of 24 and 23 cm for 2 years. Maize was not sown in 
the bare soil. This planting scheme had an actual density of 76,300 
and 82,500 plants ha~' for 2009 and 2010, respectively. Actual 
fractions of ground-mulching were ~0.5 and 0.6, respectively, 



which were defined as one minus the ratio of the summed surface 
areas of bare soil and holes to ground area. For the 0-1.0 m soil 
depth, the soil type was silt loam, with a bulk density of 1.38 g 
cm '^j a field capacity of 0.30 m"* m ^, and a wilting point of 0.12 
m'' m Over the entire growing season, maize was border- 
irrigated four times, with a total irrigation water amount of 
420 mm for both years. The amount of each irrigation event was 
measured by a water meter. Each irrigation amount was 105 mm 
on June 15, July 6, July 29 and August 20, 2009, and 105, 120, 90, 
and 105 mm on June 22, July 27, August 5 and August 29, 2010, 
respectively. 

2 Measurements of Evapotranspiration and Soil 
Evaporation 

XET was measured using an eddy covariance (EC) system 
installed in the center of the maize field. The EC consists of a fast 
response 3D sonic anemometer (CSAT3, Campbell Scientific Inc., 
UT, USA), a Krypton hygrometer (KH20, Campbell Scientific 
Inc.) and a temperature and humidity sensor (HMP45C, Vaisala 
Inc., Helsinki, Finland). All sensors were connected to a data 
logger (CR5000, Campbell Scientific, Inc.). The sonic anemom- 
eter and Krypton hygrometer were installed at a height of 1.0 m 
over the crop canopy. Net radiation (i?„) was measured by a net 
radiometer (NR-LITE, Kipp & Zonen, Delft, Netherlands), 
installed at a height of 3.5 m. Two soil heat fluxes (HFPOl, 
Hukseflux, Delft, Netherlands) were installed at a soil depth of 
8.0 cm under the plastic film and bare soil. Soil temperature above 
each soil heat flux plate was measured using thermocouples at 
depths of 0.0 cm, 2.0 cm and 6.0 cm. Soil water content from 0- 
10.0 cm was measured using a soil moisture reflectometer 
(EnviroSMART, Sentek Sensor Technologies, SA, Australia), 
Ground heat flux (G) was estimated by correcting heat fluxes at 
8.0 cm for heat storage above the transducers. The heat storage 
was determined from changes in soil temperature and moisture 
above the transducers. Based on the covariance of the 10 Hz air 
temperature and specific humidity with vertical wind velocit)', the 
latent heat flux in 30 min durations was computed using the eddy 
covariance methodolog}' with the CarboEurope recommcuidations 
[36]. Daytime XETwas adjusted by the Bowen-ratio forced closure 
method, and nighttime aET was adjusted using the filtering 
interpolation method as proposed by Ding et al. [37]. 

Soil evaporation {E^ was measured by the micro-lysimeter. 
Eight micro-lysimeter cylinders, made from PVC tubes with a 
diameter of 10 cm and height of 20 cm, were installed in bare soil 
between two plastic film rows. The cylinders were weighed at 
20:00 every day by an electric scale with a precision of 0.1 g. The 
micro-lysimeters were reinstalled within one day after each 
irrigation and hea\'y rain. E, at the field scale can be calculated 
by weighting the fraction of ground-mulching from the 
following equation. 



AM, 



(46) 



where AMj is the mean weight change of micro-lysimeter every 

day. A,, is the cross sectional area of the micro-lysimeters (78.5 cm^ 
here), is water density (1.0 g cm ''') and 10 is a conversion 
factor for changing units from cm to mm. 

3 Other IVleasurements 

Solar radiation, precipitation, air temperature, relative humidity 
and wind speed were measured with a standard automatic weather 
station (Hobo, Onset Computer Corp., USA) at a height of 2.0 m 



PLOS ONE I www.plosone.org 



6 



April 2014 | Volume 9 | Issue 4 | e95584 



Scaling Up Stomatal Conductance from Leaf to Canopy 




10:00 12:00 14:00 16:00 18:00 08:00 10:00 12:00 14:00 16:00 18:00 
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Figure 4. Diurnal variations of half-liourly canopy stomatal conductance (G^^) calculated by the big-leaf IG^ci) and dual-leaf models 
{Gsc2), and measured G^c {Gscm) inverted by the S-W model, respectively for four typical clear-sky days in different growth stages of 
maize in 2009. Leaf area index (MO was 0.16, 2.62, 5.38 and 2.99, on (a) May 22, (b) June 23, (c) July 27, and (d) September 17, respectively. 
doi:l 0.1 371/journal.pone.0095584.g004 



above the ground. Volumetric soil water content in the root zone 
was measured with PVC access tubes using the portable 
device Diviner 2000 (Sentek Sensor Technologies). Measurements 
were made at intervals of 0.1 m with a maximal soil depth of 
1.0 m at intervals of 3-5 days. Additional samplings were 
conducted before and after irrigation events, as well as after 
rainfall events. The measurements were calibrated by oven drying 
of sou samples. Interpolation was apphed between consecutive 
irrigations to determine 9^ for each day. In addition, two sets of 
ECH20 probes (Decagon Devices Inc., Pullman, WA, USA) were 
added to monitor .soil moisture at 30 min intervals in 2010. 

Ten maize plants were randomly selected to measure leaf length 
and width, and height at intervals of ~ 1 0 days during the growing 
season. Leaf area was calculated by summing the rectangular area 
of each leaf (leaf length x maximum width) multiplied by a factor of 
0.74, a conversion factor obtained by analyzing the ratio of the 
rectangular area to the real area measured by an AM300 (ADC 
BioScientific Ltd., UK). LAI is defined as maize green leaf area per 
unit ground area. The daily LAI was obtained by linear 
interpolation. 

Leaf-scale physiological measurements were performed to 
derive the stomatal conductance model parameters. A LI-6400 
portable photosynthesis system (Li-Cor Inc., Lincoln, NE, USA) 
was used to measure leaf stomatal conductance on the first fuUy 
expanded leaf, which is the fourth leaf counted from the top of the 
shoot. The diurnal measurements of leaf gas exchange were 
performed once every 2 h from 8:00 to 18:00 on six sunny days in 
10 randomly sc'l(;cted maize plants. Care was taken to keep leaves 
in their natural positions during measurement. The response of 
leaf stomatal conductance to varying PAR was measured at 30°C 
and at a CO2 concentration of 400 mol mol ^ on 29 August, 
2009. Measurements were taken at PAR levels of 2000, 1600, 
1300, 1000, 800, 600, 400. 200, 100, ,50, 20, and 0 |xmol m"^ s"'. 
The stomatal light-response curve was fit by a rectangular 



h^'perbtila to obtain the parameter values of kQ_ using the Jarvis- 
Stewart model. 

4 Model Performance 

Half-hourly G„ and XETwere calculated using the big-leaf and 
dual-leaf models Eqs. (39) and (40) with the dual-source equation 
based on the half-hourly measured meteorological data. The LAI 
and soil water were set as constants at the half-hourly time scale. 
Daily XETwas calculated using Eqs. (39) and (40) with the dual- 
source equation based on the measured average daily meteoro- 
logical data. We evaluat(;d the two models by comparing with 
measurements taken over an irrigated maize field. 

The parameters in the Jarvis-Stewart model were obtained 
using measurements of the stomatal light-response curve and the 
diurnal leaf gas exchange calculated by non-linear least-squares 
analysis (SPSS 13.0, SPSS Inc., Chicago, IL, USA). There were 
~15 days with no crop cover before the emergence of maize, 
providing the opportunity to parameterize the empirical coeHi- 
cients in the soil surface resistance model using the flux 
observations. 

The coefficient of determination (R^), root mean square error 
(RMSE) and the WiUmott's index of agreement (d) were used to 
evaluate model performance [38]. 

Results 

1 Model Parameter Estimation and Sensitivity 

From the stomatal light-response curve, we obtained best-fitting 
estimates of A:Q^by non-linear least-squares analysis using the Jarvis- 
Stewart model (Table 1). The stress coefficients of VPD and 9e, kn 
and k^, were optimized using the diurnal measurements of leaf gas 
exchange (Table 1). Soil surface resistance (r,,) was calculated by 
inverting the flux-resistance equation for the case of no crops [12]. 
Based on the relationship between ?■„ and relative soil water 
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content {6^/6,} of the top soil, the best-fitting parameters oi bj, 

and b-j were obtained (Table 1). 

The big-leaf model estimates did not closely match the response 

of G„ to incident PAR above the canopy (Q„) predicted by the 

dual-leaf model (Fig. 1). At lower LAI (2.0), the big-leaf model 
^ overestimated G„ by up to 47.4% at an intermediate irradiance 

■I (300 W m~^). At higher L4/(5.0), the big-leaf model overestimat- 

Q. ed G,c when Q„<200 W m ^ and underestimated G„ when Q^> 

^_ 200 W m ^. The sensitivity of the dual-leaf model was further 

1^3 analyzed by investigating the variations of the sunlit and shaded 

E leaf area index {LAI,i and LAI,i^ against the different LAI (Fig. 2). 

-D LAI.d approached a maximum of 1.6 when Z/1/&3.0, while LAI,h 

^ almost linearly increased as LAI increased. As a result, the ratios of 

^ LAId and LAI,it to LAI (A,i and A,i^ nonlinearly decreased and 

5 increased, respectively, as L4/ increased. 

I The diurnal variations of modeled irradiance absorbed by the 

S entire canopy [QJ) and its separation into sunlit and shaded 

> fractions {Qjt and Q,;,) are shown in Fig. 3a at higher LAI= 5.0, 

= using the measured diurnal courses of meteorological and 

^ environmental variables from June 23, 2009. Q^, Q^i and Q^/, 

™ exhibited the typical diurnal patterns, while Qj/, had a lower 

1 magnitude than Q^/ throughout the day. The average Q^i 

^ accounted for 84.2% of the Q^. The partitioning of leaves into 

S sunlit and shaded fractions continually changes throughout the 

° day (Fig. 3b). LAI^i was a convex parabola, whUe LAI,/, was a 

I concave parabola. The magnitude of LAI,/, was greater than LAI,i 

_j even during midday when the solar zenith angle is lowest at the 

higher LAL, which is consistent with the result in Fig. 2. The 
^ diurnal distributions of G„ calculated by the big-leaf model {Gjci) 

Sj and dual-leaf model {G„2) and its separation between sunlit and 

shaded fractions {G^i and G,/,) are presented in Fig. 3c. G„/ was 
° significandy lower than G^^ through most of the day at the higher 

■S LAI. G,/, showed a pattern similar to G,;, while the magnitude of G,/, 

Z, was lower than G,i. The maximum difference between them 

I occurred at the midday around 11:00. 

I 2 Comparisons of Canopy Stomatal Conductance by Big- 

.2 leaf and Dual-leaf Models 

^ Diurnal patterns of G„ calculated by the big-leaf and dual-leaf 

models (G^^ and G„2, respectively) and the measured G„ (G^J 
° calculated by using Eq. (17) are shown in Fig. 4 for four typical 

o clear-sky days in different growth stages of maize in 2009. The 

g. diurnal courses of the G,cj and G„2 were similar to the G,j„, but 

c there were differences in magnitude. G„i was significantly lower 

^ than G„m, while G„2 closely matched G^^m at higher LAI (Fig. 4c). 

S G„i was higher than G„„ at the lower LAI where there was stUl 

^ good agreement between G„2 and G,cm values (Fig. 4b and d). 

^ Although the slope of the linear regression between G„; and G„„ 

m was 0.97, G„i was overestimated at the lower G„„ and 

§ underestimated at the higher G„m (Fig. 5a), as shown by the lower 

c r2 of 0.81, greater RMSE of 1.7047 mm s"' and lower d of 

S ^ 0.9563, indicating that the big-leaf model yielded large errors for 

E S estimating G„. In contrast, the slope of the linear regression 

^ M between G,c2 and G„m was 1.01, with an R^ of 0.98, RMSE of 

S '& 0.6120 mm s~' and d of 0.9951, indicating that there was good 

data-model agreement between predictions and measurements 
(Fig. 5b). 

3 Comparisons of Crop Evapotranspiration by Big-leaf 
and Dual-leaf Models 

Diurnal variations of half-hourly kET calculated by the dual- 
source big-leaf [XET]) and dual-leaf models {XET2}, and measured 
aET [aET^, respectively are presented for four typical clear-sky 
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Figure 5. Relationship between canopy stomatal conductance [G^^ estimated by the big-leaf [Gsc) and dual-leaf models {Gscj), and 
measured G^c [Gscm) Inverted by the S-W model, respectively for four typical clear-sky days in different growth stages of maize in 
2009. 

doi:1 0.1 371 /journal.pone.0095584.g005 



days in different growth stages of maize in 2009 (Fig. 6). The 
diurnal patterns of estimated XET were similar to the measure- 
ments. aETj wa.s overestimated at lower LAI, and underestimated 
at higher LAI. The linear regression presented that aETj was 
overestimated by 8.7% (R^ = 0.97) and 19.7% (R^ = 0.96) for 
Z^/=2.62 and 2.99, respectively (Fig. 6b and d). lETj was 
underestimated by 13.9% (R^ = 0.97) for LAI=538 (Fig. 6c). In 



contrast, XET2 had a good agreement with measurements for 
differing LAI, with a linear slope of 1.01 and an of 0.97. 

The irrigation scheduling and mulching fractions in 2009 and 
2010 were different (See section 3.1), which yielded different LAI 
and soil water regimes for the two years (See Figure 1 and Table 1 
in Ding et al.[39]). The maximum and averaged values of LAI 
respectively were 5.4 and 3.1 for 2009, and 4.7 and 2.7 for 2010. 
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Figure 6. Diurnal variations of half-hourly crop evapotranspiration (AET) calculated by the dual-source with the big-leaf (ifr,) and 
dual-leaf models (AET^, and measured 2ET[2ETni, respectively for four typical clear-sky days in different growth stages of maize in 

2009. Leaf area index {LAI) was 0.16, 2.62, 5.38 and 2.99, on (a) May 22, (b) June 23, (c) July 27, and (d) September 17, respectively. 
doi:1 0.1 371/journal.pone.0095584.g006 
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Figure 7. Comparison between Kialf-Kiourly estimated evapotranspiration by tlie dual-source witli tlie big-leaf UET,) and dual-leaf 
models UET^ versus measured UETby eddy covariance (JIET^ during the entire growth period of maize in (a and b) 2009 and (c and 
d) 2010. 

doi:! 0.1 371/journal.pone.0095584.g007 



The extractable soil water in the root zone [Be) was significantly 
different between the two years. Before the first and second 
irrigation events in 2010, there were 9 and 12 days of Be below 
50% total available water in the root zone (TAIV), which was 
regarded as a threshold of crop water stress [13]. Conversely, most 
Be WHS higher than 50% of TAW in 2009. All of these differences 
led to differing 2.ET and its components, which provided a good 
dataset to test the big-leaf and dual-leaf models over two different 
hydrometeorological and management strategies. 

The scatterplots of half-hourly XET exhibited that aETj was 
overestimated for lower values and underestimated for higher 
values, respectively (Fig. 7a and 7c). Total ^ETj was underesti- 
mated, with a slope of linear regression of 0.94 (S3 = 0.83) and 
0.93 (R^ = 0.82), respectively, for 2009 and 2010. RMSE was 
72.22 and 70.97 W m"^ and d was 0.9521 and 0.9472 for 2009 
and 2010, respectively (Table 2). In contrast, there was good data- 
model agreement between measurements and estimated half- 
hourly AETz in 2009 and 2010 (Fig. 7b and 7d). The slopes of 
linear regressions between the estimates and measurements were 
1.02 and 1.03, with of 0.90 and 0.88, RMSE of 58.06 and 
62.31 W m"2 and d of 0.9706 and 0.9626 for 2009 and 2010, 
respectively (Table 2). Daily estimated aET^ enhanced data-model 
agreement, with R'^ of 0.91 for the 2 years despite the linear slopes 
were the same as those of half-hourly values (data not shown). The 
statistical test showed that the slopes were not significanfly 
different with one (P=0.114 and 0.092), and the intercepts were 
not significantly different with zero (P = 0.2 1 5 and 0. 1 74) for the 2 
years. 



Seasonal variations of daily estimated and measured E^ using 
Eq. (15) combined with Eqs. (39) and (40) are presented in Fig. 8 
for 2009 and 2010. Both the dual-source big-leaf and dual-leaf 
models could capture the variability of E, even when irrigation or 
precipitation occurred and when the canopy partially covered the 
ground during the initial growth stages. However, daily values of 
E^] were underestimated at the early and late stages and 
overestimated at the middle stage (Fig. 8a and b). In general, 
total E,j was underestimated by 5% and 7% for 2009 and 2010, 
respectively (Fig. 8c and d). In contrast, there was satisfactory data- 
model agreement between predicted and measured E^ using the 
DSDL model for the two years. The slopes of linear regressions 
between estimates and measurements were 1.02 (R2 = 0.68) and 
1.01 (R2 = 0.73), with RMSE of 0. 1220 and 0.1565 mm d"' and d 
of 0.9129 and 0.9218 for 2009 and 2010, respectively (Table 2). 

Discussion 

In this paper, we have extended the big-leaf model by 
developing a dual-leaf model. The dual-leaf model presented here 
is an improvement over the previous big-leaf model, as more 
realistic non-uniform vertical profiles of radiation and stomatal 
conductance are now incorporated into the model. The penetra- 
tion of beam radiation, its variation and the dynamics of sunlit and 
shaded LAI throughout the day all affect the ability of the big-leaf 
model to simulate diurnal changes in Gsc [18,29]. In the dual-leaf 
model, these canopy features can be expUcitly incorporated by 
dividing the canopy into sunlit and shaded fractions and modeling 
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Figure 8. Comparison between daily estimated soil evaporation (E^^ by the dual-source with the big-leaf (f^,) and dual-leaf models 
(Esil versus measurements (£5^) during the entire growth period of maize in 2009 and 2010. Seasonal variations of the estimated and 
measured fj against days after sowing (DAS) are presented in (a) and (b). The linear regressions between them are presented in (c) and (d). 
dol:l 0.1 371/journal.pone.0095584.g008 



each fraction of G^^ by scaling up the respective stomatal 
conductance separately. It is more complex than the big-leaf 
model, but the dynamic partitioning of LAI and irradiance 
between sunlit and shaded leaves has further reduc(-d the errors 
associated with simplifying the leaves to only a big-leaf using either 
the total or empirically effective LAI. 

The ability of the dual-leaf model was examined by comparing 
the estimated values and actual measurements. In contrast to 
significant errors by the big-leaf model, the dual-leaf model 
accurately reproduced the variation of G^c (Fig. 4 and 5). One 
reason the dual-leaf model works so well is that it accommodates 
the nonlinear response of stomata to light [25,40] . Stomata-light 
responses of leaves can vary with depth in the canopy and this 
variation can be incorporated by partitioning the canopy into 
several layers and estimating the sunlit and shaded leaf fractions in 
each layer [20,34]. However, usually this is not necessary and a 
single, representative light response curve can be used for the 
entire canopy [17,40]. In the dual-leaf model, the entire G„. may 
be calculated by summing contributions of sunlit and shaded 
leaves. These two contributions are added separately because 
sunlit leaves wUl be light-saturated while shaded leaves will be in 
the linear portion of the light-stomata relationship [18,33]; thus Gsc 
is not proportional to average light levels [18]. Because of the 
nonlinear relationship between stomatal conductance and PAR, 
the predicted G^^ will be overestimated when the average absorbed 
PAR is used to scale up the leaf stomatal conductance for the 
entire canopy as in the big-leaf model [17,40]. On the other hand, 
the G^c is underestimated when the effective LAI h used to scale up 
the leaf stomatal conductance at higher LAI ((Fig. 4c). Since the 



nonlinear relationship was considered and the sunlit-shaded 
method was introduced in the dual-leaf model, estimates using 
the dual-leaf model closely match the measured G„, (Fig. 4 and 5). 

The performance of the DSDL for estimating ).ET was 
investigated. The DSDL is robust for estimating XET over a 
range of canopy leaf areas and environmental variables (Fig. 6, 
Fig. 7 and Table 2). The good data-model agreement indicated the 
strengths of the DSDL model as a model framework and the 
reference for validating other approaches of calculating the G,£ 
using the dual-leaf model. Our framework of modeling XET also 
provided a soil evaporation estimation model Eq.(15) with a 
modified soil surface resistance term, which is useful for enhancing 
crop production by reducing the E^ fraction of kET [10,15]. 

The DSDL model is physically process-based, yet sufficiently 
simple to be effectively parameterized. The dual-leaf model 
requires only four additional equations, Eqs (32), (33), (43) and 
(44), beyond those required in the model of leaf stomatal 
conductance, to calculate the LAI and absorbed irradiance of the 
sunlit and shaded leaves. This simplicity makes it attractive for 
incorporation into in crop models, land surface schemes, and 
regional or global water cycle studies [29,40] . This model can also 
be used to assess effects of climate change on crop ecophysiology. 

Conclusions 

In this paper, a dual-leaf model for scaling-up stomatal 
conductance from the leaf to the canopy level was developed 
through the dynamic partitioning of the leaf area index and 
irradiance between sunlit and shaded leaves. In the model, canopy 
stomatal conductance was calctilated by dividing the canopy into 
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sunlit and shaded fractions and each fraction was modeled 
separately based on the absorbed irradiances. The dual-leaf model 
provided estimates of G^^ which were nearly the same as 
measurements, and were significantly more accurate than those 
of the big-leaf model. Our results showed excellent agreements 
between XET measurements gathered by the eddy covariance 
technique over an irrigated maize field during 2009 and 2010 
under two different hydrometeorological and management condi- 
tions, and estimates of lET using the DSDL. The framework of 
the model can also satisfactorily estimate soil evaporation. Our 
proposed model provides an alternative approach to calculate 
lET^ which is simple and attractive for incorporation into other 
comprehensive crop models. 
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